function err_mom = objective(x, Markuptarget)



printr       = 1; 

if printr == 1
    
    options   = optimoptions('fsolve', 'FiniteDifferenceStepSize', 1e-4, 'Display', 'iter', 'TolX', 1e-11, 'TolFun', 1e-11, 'UseParallel', false);

else
    
    options   = optimoptions('fsolve', 'FiniteDifferenceStepSize', 1e-4, 'Display', 'none', 'TolX', 1e-11, 'TolFun', 1e-11, 'UseParallel', false);

end


% Calibrated Parameters

N            = x(1);                         % mean number of firms in sector
p.sigma      = x(2);
p.gamma      = x(3); 
p.xi         = x(4); 


% Assigned Parameters

p.beta       = 0.96;                        % period is 1 year
p.varphi     = 0.04;                        % exit rate 
p.delta      = 0.06;                        % depreciation rate

p.r          = 1/p.beta - 1; 
p.nu         = 1;                           % inverse labor supply elasticity
p.alpha      = 1/3;                         % capital share in value added
p.theta      = 0.5;                         % elasticity of substitution between VA and intermediates
 
S            = 25000;                       % number of sectors
Nmax         = max(100, floor(N*2));        % maximum number of firms in sector



% Pareto tail productivity

Y            = 1; 
R            = p.r + p.delta; 

mshare       = 0.45;                        % share of intermediate inputs in total sales


% Quadrature
 
 
nz          = 100;                          % Gaussian discretization for z

zmax        = -1/p.xi*log(1e-12);


[p.z, p.w]  = qnwlege(nz, 0, zmax); 

p.w         = p.xi.*p.w.*exp(-p.xi.*p.z);

p.w         = p.w./sum(p.w); 
p.z         = exp(p.z); 

cut         = cumsum(p.w) > 1 - 1e-5; 

wend        = sum(p.w(cut)); 
zend        = p.w(cut)'*p.z(cut)/wend; 

p.w         = [p.w(~cut); wend]; p.w = p.w/sum(p.w);
p.z         = [p.z(~cut); zend]; 

nz          = size(p.z, 1); 
Fcum        = [0; cumsum(p.w)];              % cumulative distribution

Fcum(end)   = 1; 

% Random draws

rng(0); 

% 1. Allocate firms to sectors

ns          = max(1, poissrnd(N, S, 1)); 
Nf          = Nmax*S; 

% 2. Productivity draw from discretized distribution

unif        = nodeunif(Nf, 0, 1); 
unif        = unif(randperm(Nf)); 

[~, ~, bin] = histcounts(unif, Fcum);

z           = reshape(p.z(bin), Nmax, S); 
zsave       = z; 
nssave      = ns; 

Nsave       = mean(ns); 




for i = 1 : S
   
    z(ns(i) + 1 : end, i) = eps;         % kill extra firms (set z = eps)
    
end

nsmax         = max(ns); 
z             = z(1 : nsmax, :);         % no reason to carry extra entries

z             = sort(z, 1, 'descend');


% Solve equilibrium markups in every sector

[mu, prof, omega, Mus, Zs, Omegas, Mu, Z, Omega]  = industryequilibrium([], z, p);

p.phi            = 1 - mshare*Omega^(1 - p.theta)*Mu;                                       % choose to match materials share in sales = mshare 

W                = (1 - p.alpha)*(((Omega^(1 - p.theta) - (1 - p.phi))/p.phi)^(1/(1 - p.theta))*(R/p.alpha)^(-p.alpha))^(1/(1 - p.alpha));        % Invert the CES Cost aggregator                                              % Labor only factor of production                

Lp               = (1 - p.alpha)*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/W; 
K                =       p.alpha*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/R;
B                = (1 - p.phi)/Mu*(1/Omega)^(1 - p.theta)*Y;

C                = Y - p.delta*K - B;


moment_data        = [Markuptarget; 0.45; 0.430; 0.723; 0.213]; 


moment_model       = zeros(size(moment_data)); 

moment_model(1)    = Mu; 
moment_model(2)    = B/Y;


sales              = omega.*Omegas*Y;
labor              = (1 - p.alpha)*p.phi*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*sales./mu;

lshare             = sum(labor)'./sum(sales)';               % sectoral labor share 
hh                 = sum(omega.^2)'; 

bhat               = lscov([ones(S, 1), hh], lshare);

CR4                = Omegas*sum(omega(1:4, :))'/sum(Omegas);
CR20               = Omegas*sum(omega(1:20, :))'/sum(Omegas);

moment_model(3)    = CR4;
moment_model(4)    = CR20;
moment_model(5)    = -bhat(2);

if printr
    

fprintf('\n');

fprintf('Moments used in Calibration: Left = Model, Right = Data  \n');   
fprintf('\n');
fprintf('\n');
fprintf('Aggregate Markup                           = %9.3f %9.3f\n',  [moment_model(1),    moment_data(1)]);
fprintf('Share Intermediates in Sales               = %9.3f %9.3f\n',  [moment_model(2),    moment_data(2)]);
fprintf('\n')
fprintf('CR4                                        = %9.3f %9.3f\n',  [moment_model(3),    moment_data(3)]);
fprintf('CR20                                       = %9.3f %9.3f\n',  [moment_model(4),    moment_data(4)]);
fprintf('\n');
fprintf('Regression coefficient                     = %9.3f %9.3f\n',  [moment_model(5),    moment_data(5)]);
fprintf('\n');

end



weights        = ones(numel(moment_model), 1); 

weights([1]) = 100; 
weights([5]) = 10; 

weights        = weights/sum(weights);


err_mom        = (moment_model - moment_data)./(1 + moment_data);
err_mom        = (weights'*err_mom.^2).^(1/2);

if printr
        
fprintf('\n');   
fprintf('Error in Moments                           = %9.6f \n', err_mom);
fprintf('\n');
fprintf('\n');   

end

format short g
fprintf('\n');
disp([x(:)',  err_mom]);
